A novel liver zonation phenotype-associated molecular classification of hepatocellular carcinoma

Background Liver zonation is a unique phenomenon in which the liver exhibits distinct functions among hepatocytes along the radial axis of the lobule. This phenomenon can cause the sectionalized initiation of several liver diseases, including hepatocellular carcinoma (HCC). However, few studies have explored the zonation features of HCC. Methods Four single-cell RNA sequencing datasets were used to identify hepatocyte-specific zonation markers. Integrative analysis was then performed with a training RNA-seq cohort (616 HCC samples) and an external validating microarray cohort (285 HCC samples) from the International Cancer Genome Consortium, The Cancer Genome Atlas, Gene Expression Omnibus, and EMBL’s European Bioinformatics Institute for clustering using non-negative matrix factorization consensus clustering based on zonation genes. Afterward, we evaluated the prognostic value, clinical characteristics, transcriptome and mutation features, immune infiltration, and immunotherapy response of the HCC subclasses. Results A total of 94 human hepatocyte-specific zonation markers (39 central markers and 55 portal markers) were identified for the first time. Subsequently, three subgroups of HCC, namely Cluster1, Cluster2, and Cluster3 were identified. Cluster1 exhibited a non-zonational-like signature with the worst prognosis. Cluster2 was intensively associated with a central-like signature and exhibited low immune infiltration and sensitivity toward immune blockade therapy. Cluster3 was intensively correlated with a portal-like signature with the best prognosis. Finally, we identified candidate therapeutic targets and agents for Cluster1 HCC samples. Conclusion The current study established a novel HCC classification based on liver zonation signature. By classifying HCC into three clusters with non-zonational-like (Cluster1), central-like (Cluster2), and portal-like (Cluster3) features, this study provided new perspectives on the heterogeneity of HCC and shed new light on delivering precision medicine for HCC patients.


Introduction
Hepatocellular carcinoma (HCC) is one of the most common types of cancer worldwide (1). Despite the rapid progression of new diagnostic methods and therapeutic strategies for HCC, its prognosis is still unfavorable due to its high heterogeneity. Therefore, uncovering the molecular mechanisms underlying HCC diversity is essential for the development of targeted and effective therapies.
The liver is a central organ that maintains physiological homeostasis. Liver lobules are the functional units of the liver and are hexagonal in shape. Hepatocytes are organized in a cord-like arrangement along the liver lobule, extending from portal nodes (PN) to the central vein (CV) (2). According to the location of hepatocytes, the liver lobule can be divided into three zones: zone 1 (periportal area) is the region near the portal triad, zone 3 (pericentral zone) is near the CV, and the region between these zones is zone 2 (mid-lobular) (3). The liver shows functional partition along the lobule radial axis, a phenomenon known as "liver zonation". Hepatocytes in different zonations perform various functions, allowing multiple of functions to proceed in parallel. Most of the liver's metabolic functions take place in zone 1 hepatocytes, such as b-oxidation, gluconeogenesis, and lipid metabolism. On the contrary, zone 3 hepatocytes play central roles in glycolysis, xenobiotic biotransformation reactions, and glutamine synthesis (4).
Zonation is a dynamic process that plays essential roles in regulating liver disease phenotypes and progression. For example, since fatty acid synthesis and lipid accumulation occur predominantly in zone 3, whereas fatty acid oxidation occurs in zone 1. Patients with non-alcoholic fatty liver disease (NAFLD) often (approximately 37%) exibit perivenous dominant steatosis (5). In addition, patients with HCC exhibit impaired Wnt/b-catenin signaling (6), which stands out as a major regulator of liver zonation, regulating about a third of liver zonated genes (7).
In the past decades, genome-wide analyses have been devoted to deciphering the molecular mechanisms of HCC diversity (8)(9)(10)(11)(12). However, few studies have explored the molecular classification of HCC associated with liver zonation characteristics (11). In the present study, we constructed a novel molecular classification associated with liver zonation phenotype, and 3 subgroups of HCC, namely Cluster1, Cluster2, and Cluster3 were identified. We then evaluated the prognostic value, clinical characteristics, transcriptome and mutation features, immune infiltration, and immunotherapy response of the HCC subclasses. Finally, we identified candidate therapeutic targets and agents for Cluster1 HCC samples, which have the worst prognosis among the three clusters.
RNA sequencing data (counts) of 374 and 242 HCC human samples with available clinical information were obtained from the TCGA-LIHC cohort (https://xenabrowser.net/datapages/) and the LIRI-JP cohort (https://dcc.icgc.org/projects/LIRI-JP), respectively. The SVA R package was utilized to merge and remove the batch effects of the two RNA-seq datasets to create one metadata cohort using the Sangerbox online tool (13). Additional microarray data of 225 and 60 HCC samples from GSE 14520 ( https:/ / www.ncbi.nlm.nih.gov/geo/) and E-TABM-36 (https:// www.ebi.ac.uk/arrayexpress/files/E-TABM-36) based on the HG-U133A platform, were used for external validation. A metadata cohort was also created by merging the two microarray datasets, and batch effects were removed using the combat function in SVA R package. A boxplot before and after batch effect correction is shown in Supplementary figure S1.
Moreover, gene somatic mutation data (MAF files) of the TCGA-LIHC and LIRI-JP cohorts were obtained from TCGA and ICGC databases, respectively. Additionally, copy number data of GISTIC2 for the TCGA-LIHC cohort were accessed from the UCSC Xena (https://xenabrowser.net/datapages/). Furthermore, the gene expression data and drug sensitivity data (AUC values) of hepatocellular carcinoma cell lines were downloaded from the DepMap database (https://depmap.org/ portal/download/).

Identification of hepatocyte-specified liver zonation markers
The procedure for identifying hepatocyte-specified liver zonation markers was shown in Figure 1 To identify hepatocyte markers in both mice and humans, we utilized the Seurat R package (16) and applied standard downstream processing to mouse and human liver snRNA-seq data obtained from GSE192742 (17). Firstly, the LogNormalize method was applied for data normalization. Then, the "FindVariableFeatures" function (selection.method = "vst", nfeatures = 3000) and the "FindIntegrationAnchors" function were using to select features and anchors for further integration. After integration, we eliminated low-quality cells with ≤ 500 genes or a mitochondrial gene ratio ≥ 5%. Finally, 25470 human liver cells and 25335 mouse liver cells that were deemed to be of good quality underwent further analysis. Following data integration and scaling, Seurat's "RunPCA" function and "FindClusters" function were employed to perform principal component analysis and clustering. Dimensionality reduction was conducted utilizing the "RunUMAP" function, which utilized the top 30 calculated dimensions and a resolution of 0.5. Subsequently, the "FindAllMarkers" function was used to identify marker genes for each cluster (adjusted P < 0.05). Afterward, ToppGeneSuit (https:// toppgene.cchmc.org/), CellMarker (http://xteam.xbio.top/ CellMarker/), and PanglaoDB (https://panglaodb.se/search.html) databases were employed for the annotation of cell types. Finally, the "FindAllMarkers" function was employed to identify hepatocyte markers by comparing hepatocytes and other cell types (adjusted P < 0.05 and |logFC| > 0.25).

Identification of HCC subclasses
HCC samples of training and validating cohorts were classified using non-negative matrix factorization (NMF) clustering based on 94 liver zonation markers (18). And then, the prognosis of different subclasses in the training and validation cohorts was evaluated using the Kaplan-Meier log-rank test. Moreover, subclass mapping (SubMap) analysis (19) (Gene Pattern), a method for assessing the comparability of molecular classes between different patient cohorts based on their expression patterns, was then employed to verify whether the subclusters identified in the two cohorts were associated.

Molecular characteristics of HCC subclasses
The differentially expressed genes (DEGs) among the three clusters were identified using the "limma" package in R with cutoff criteria of |log2 fold change (FC)| > 1 and an adjusted P value < 0.05. Only genes that showed substantial differences in expression across all three comparisons were designated clusterspecific DEGs. Functional enrichment analysis of cluster-specific DEGs was conducted by the Metascape database (https:// metascape.org/). P<0.05 was considered statistically significant.
Gene set variation analysis (GSVA) was then employed to estimate the score of the 50 hallmark gene sets that were achieved from the MsigDB database (http://www.gsea-msigdb.org/gsea/ msigdb/index.jsp) (20). After that, we used the "ComplexHeatmap" R package to display distinct pathways among the three HCC clusters.
Moreover, Nearest template prediction (NTP) analysis (Gene Pattern modules) was used to predict the correlation between previously published HCC molecular classifications and our classification.

Estimation of immune infiltration and prediction of the immunotherapy response
To further explore the difference in immune cell infiltration among the clusters, the CIBERSORT algorithm (18) was used to estimate the fraction of 22 immune cell types in the HCC samples using Sangerbox online tool (13). In addition, single-sample GSEA (ssGSEA) was also used to estimate immune infiltration, which computed an enrichment score representing the degree to which genes in 28 immune cell gene sets were coordinately up or downregulated within a single sample (21).
To predict the response of immunotherapy in our subclasses, SubMap analysis (Gene Pattern) was employed to compare the similarity of gene expression profiles between our subclusters and a cohort of melanoma patients with programmed cell death protein-1 (PD1) inhibitor or cytotoxic T-lymphocyte-associated protein-4 (CTLA-4) inhibitor treatment (22).

Identification of potential drug targets and therapeutic agents for Cluster1 of HCC
The cluster1-specific DEGs with overexpression (log2FC > 1 and adjusted P < 0.05) in HCC tissues and low CERES scores (< -0.5) were defined as potential drug targets. The CERES scores were acquired from the dependency map (DepMap) portal (https:// depmap.org/portal/). The CERES score is utilized to evaluate the dependency of the interest gene in a certain cancer cell lines, and a lower score suggests a higher likelihood that the gene is crucial for cell growth and survival of a given cancer cell line.
The Genomics of Drug Sensitivity in Cancer (GDSC) (23), and PRISM (24) databases contain information on drug sensitivity and gene expression profiles of cancer cells, which can be employed to establish a prediction model of drug response. The oncoPredict R package is used for predicting drug response through a ridge regression model to calculate sensitivity scores (low sensitivity score indicates high drug sensitivity) (25).

Statistical analysis
All computational and statistical analyses were performed using R programming (https://www.r-project.org/) and SPSS 22.0 (IBM Corp., Armonk, NY, USA). To compare two or three groups with normally distributed variables, the unpaired Student's t-test and one-way ANOVA were utilized, respectively. Kruskal-Wallis tests were employed to compare three groups with non-normal distribution parameters. The chi-square test or Fisher's exact tests were used for analyzing contingency table variables. Survival analysis was performed using Kaplan-Meier methods with the log-rank test. A two-tailed P value < 0.05 was statistically significant.

Identification of hepatocyte-specified zonation markers
A flowchart was depicted to systematically describe the processes of identifying hepatocyte-specified zonation markers ( Figure 1A). Firstly, we identified 1467 liver zonation-associated genes that were overexpressed in the CV or PN area in Halpern et al.'s cohort (Table S1). Moreover, 1034 liver zonation-associated genes were screened in Moshe et al.'s cohort (Table S2) (Tables S3, S4). Afterward, 105 mouse hepatocyte-specified zonation markers were identified ( Figure 2C). Finally, after transferring the 105 mouse zonation markers to human homologous genes, we screened 94 human hepatocyte-specified zonation markers, including 39 central markers and 55 portal markers ( Figure 2D, Table S5). The expression pattern of the hepatocyte-specified zonation markers in the liver lobule was shown in Figures  Consistent with previous studies, functional enrichment analysis of these zonation markers indicated that central markers were mainly significantly enriched in processes of fatty acid metabolism, bile secretion, cholestasis, xenobiotic transport, glucose metabolism, and lipid localization. However, portal markers were mainly significantly enriched in processes of amino acid metabolism, complement and coagulation cascades, the urea cycle, regulation of proteolysis, gluconeogenesis, and the response to metal ion (Table S6).

NMF identifies three clusters in HCC
The procedure of our bioinformatics analyses is shown in Figure 1B. The RNA-seq cohort comprising 616 HCC samples from TCGA-LIHC and ICGC-LIRI-JP were clustered based on the expression profile of 94 liver zonation markers using NMF consensus clustering. As shown in Figure 3A, three distinct clusters were identified in the RNA-seq cohort: Cluster 1 with 123 cases, Cluster 2 with 244 cases, and Cluster 3 with 249 cases. Moreover, a significant prognostic difference was observed among the three clusters in the RNA-seq cohort. Patients in Cluster1 had the worst prognosis, while patients in Cluster3 had the best prognosis (Overall log-rank test P = 1.041×10-6, Cluster1 vs 2 P = 1.189×10-2, Cluster1 vs 3 P = 4.275×10-6, Cluster2 vs 3 P = 4.234×10-3; Figure 3B). Subsequently, we explored the differences in the expression patterns of liver zonation markers among these clusters. As shown in Figure 3C, Cluster3 had the highest expression level of portal markers, while Cluster1 had the lowest expression level of both central and portal markers. Moreover, we performed GSVA analysis of gene sets based on central and portal markers to depict the zonation characteristics of each cluster. As shown in Figure 3D, consistent with gene expression level, Cluster3 had the highest GSVA score for portal markers, while Cluster1 had the lowest GSVA scores for both central and portal markers.
Subsequently, we conducted another independent analysis on a microarray cohort with 285 HCC samples from GSE14520 and E-TABM-36, the results of which also showed that there were three distinct clusters of HCC: Cluster 1 with 85 cases, Cluster 2 with 106 cases, and Cluster 3 with 94 cases (Supplementary figure S2A). Moreover, consistent with the RNA-seq cohort, a significant prognostic difference was observed among the three clusters in the microarray cohort, which patients in Cluster1 had the worst prognosis, and patients in Cluster2 and Cluster3 had better prognoses (Overall Survival: Overall log-rank test P = 1.394×10-4, Cluster1 vs 2 P = 1.394×10-4, Cluster1 vs 3 P = 7.750×10-4, Cluster2 vs 3 P = 5.242×10-1; Figure 3E; Recurrence Free Survival: Overall log-rank test P = 3.616×10-4, Cluster1 vs 2 P = 2.201×10-3, Cluster1 vs 3 P = 3.367×10-3, Cluster2 vs 3 P = 9.369×10-1; Figure 3F). Furthermore, as shown in Supplementary figure S2C, consistent with the RNA-seq cohort, Cluster3 in the microarray cohort had the highest expression level and GSVA score for portal markers, while Cluster1 had the lowest expression level and GSVA scores for both central and portal markers. Finally, a SubMap analysis was performed to determine whether the clusters identified in the two above datasets were associated, and the result showed that Cluster1, Cluster2, and Cluster3 in the RNAseq cohort were highly associated with corresponding clusters in the microarray cohort (P = 0.009), indicating there were three distinct molecular subclasses of HCC with different gene expression patterns (Supplementary figure S2D).

Correlation of the zonation-associated clusters with clinical characteristics and molecular subclasses of HCC published previously
We then explored the association between HCC-related clinicopathological variables and our classification based on the RNA-seq ( Figure 4A and Table S7) and microarray cohorts (Supplementary figure S3A and Table S8). The chi-square and ANOVA tests indicated significant relationships between clinicopathological characteristics and HCC subtypes in the RNAseq cohort. Favorable survival status (P < 0.001), longer survival time (P < 0.001), lower serum AFP level (P < 0.001), bigger BMI (P=0.027), lack of vascular invasion (P = 0.004), early TNM stage (P < 0.001), early histologic grade (P < 0.001), were associated with Cluster3, while worse survival status, shorter survival time (P < 0.001), smaller BMI, presence of vascular invasion, advanced histologic grade, advanced TNM stage, and high serum AFP level were associated with Cluster1. Similarly, in the microarray cohort, Cluster3 was correlated with favorable survival and recurrence status (P = 0.009, P = 0.045), longer survival and recurrence time (P < 0.001, P < 0.001), lower serum AFP level (P < 0.001), smaller tumor size (P = 0.01), and early TNM and BCLC staging (P = 0.005, P = 0.041). However, worse survival and recurrence status, shorter survival and recurrence time, smaller BMI, presence of vascular invasion, advanced histologic grade, advanced TNM stage, and high serum AFP level were associated with Cluster1.

Transcriptomes of the zonationassociated HCC clusters
Differential analyses were carried out to comprehend the distinctions in the molecular and biological processes among the three HCC subclasses. Significant differences in gene expression were defined as |log2FC| > 1 and the adjusted P-value < 0.05. Only genes that showed substantial differences in expression across all three comparisons were designated cluster-specific DEGs. Finally, 2682 specific DEGs (1868 upregulated and 814 downregulated) for Cluster1, 993 specific DEGs (118 upregulated and 875 downregulated) for Cluster2, and 405 specific DEGs (321 upregulated and 84 downregulated) for Cluster3 were identified (Table S9).
Next, functional enrichment analysis of the cluster-specific DEGs was performed utilizing the Metascape database, and remarkably enriched biological processes or pathways are shown in Table S9. The specific DEGs of the three clusters showed distinct enrichment of biological processes. Cluster1 was enriched in some differentiation and development-relevant processes. However, multiple metabolism-related biological processes and pathways were significantly enriched for Cluster3. For Cluster2, it was enriched in the processes associated with molecular transport and the WNT signaling pathway (Table S10).

The gene mutation profile of the zonation-associated HCC clusters
To further investigate the different patterns of gene mutations among HCC clusters, somatic mutation data from the RNA-seq cohort was analyzed. The top 25 genes with the highest mutation rates in all three clusters were shown in Figures 6A-C. In Cluster1, TP53 exhibited the highest somatic mutation rate (42%), followed by TTN (27%), MUC16 (17%), PCLO (14%), RYR2 (13%), and LRP1B (12%). In Cluster2, CTNNB1 exhibited the highest somatic A previous study identified 161 HCC driver genes that may play critical roles in the treatment of HCC (27). Therefore, we further explored the mutation profile of these driver genes. The top 20 driver genes with the highest mutation rates in all three clusters were shown in Figure 6D. TP53 exhibited the highest somatic mutation rate (29%), followed by CTNNB1 (25%), ALB (12%), APOB (10%), ARID1A (8%), ARID2 (8%), and AXIN1 (6%). Moreover, we compared the tumor mutation burden (TMB) of HCC driver genes among the three clusters. The results showed that Cluster2 exhibited a higher TMB level than the other two clusters ( Figure 6E).

Correlation of the zonation-associated HCC clusters with immune infiltration and immunotherapy response
Given the significant differences in immune processes among clusters, immune infiltration was examined to characterize their immunological landscape. As shown in Figure 7A, the abundance of immune cell types was calculated using CIBERSORT and ssGSEA algorithms. As shown in Figure 7B, compared with the other two subclasses, Cluster2 presented the lowest abundance of 24 immune cell types, including B cells memory, macrophages M2, neutrophils, activated B cells, central memory CD4 T cell, central memory CD8 T cell, effector memory CD4 T cell, effector memory CD8 T cell, immature B cell, T follicular helper cell, regulatory T cell, Type 1/2/ 17 helper cell, activated dendritic cell, CD56dim natural killer cell, immature natural killer cell, macrophage, mast cell, MDSC, monocyte, natural killer T cell, natural killer cell, and Plasmacytoid dendritic cell. We further investigated the differences in the expression of 47 immune checkpoint genes among the three subclasses. Consistent with immune infiltration, Cluster2 exhibited the lowest expression for 44 immune checkpoint genes (except for ICOS, NRP1, and TNFSF14) compared to Cluster1 and Cluster3 ( Figure 7C). Different immune infiltration and immune checkpoint gene expression patterns among HCC clusters indicate that the immunotherapy response among the clusters may be different. Therefore, we compared the expression profiles of the three clusters (Cluster1, Cluster2, and Cluster3) with another published melanoma cohort containing 47 patients who received PD1 inhibitor or CTLA-4 inhibitor using SubMap analysis. As shown in Figure 7D, the expression profile of Cluster2 was correlated with the CTLA4 response group (P = 0.05), while Cluster3 was associated with the PD1 response group (P = 0.003), suggesting that patients within the Cluster1 group were more likely to respond to the anti- CTLA4 therapy and Cluster3 group was more likely to respond to anti-PD1 therapy.

Identification of potential drug targets and candidate therapeutic agents for Cluster1
Since patients in Cluster1 have the worst prognosis, we further identified the potential drug targets for Cluster1. The cluster1specific DEGs with overexpression (log2FC > 1 and adjusted P < 0.05) in HCC tissues and low CERES score (< -0.5) were defined as potential drug targets. Eventually, a total of 34 potential drug targets were identified ( Figure 8A, Table S11). Among them, drugs (compounds) that target AURKB, BIRC5, KIF11, PLK1, PLK4, RAD51, TOP2A, TTK, and TUBB3 were identified ( Figure 8A, Table S11).
To further identify the potential therapeutic agents for Cluster1, the oncoPredict R package was used to calculate the sensitivity score of the drugs (compounds) for potential drug targets of Cluster1 based on the drug sensitivity and gene expression profiles of HCC cell lines from GDSC and PRISM databases. As shown in Figures 8B, C, eleven PRISM-derived compounds (including Teniposide, K-858, TAK-901, Dexrazoxane, Mitoxantrone, AT-9283, Axitinib, Podophyllotoxin, Amonafide, Amsacrine, and Idarubicin) and six GDSC-derived compounds (including Axitinib, Doxorubicin, Etoposide, GSK1070916, YM−155, and Alisertib) were found to be more sensitive in Cluster1 HCC patients.

Discussion
Liver zonation is a unique phenomenon that exhibits a distinct division of functions among hepatocytes along the lobule radial axis, optimizing overall liver function, and playing a central role in regulating liver disease phenotypes and progression, including HCC (28). Therefore, obtaining a better understanding of the characteristics of liver zonation in HCC may help us unveil the mechanisms of HCC heterogeneity and develop more effective therapeutics for HCC. With the breakthrough of RNAsequencing technology, scRNA-seq coupled with spatial mapping has demonstrated previously unknown molecular patterns of hepatocytes, shedding new light on the functional features of hepatocytes across different zones of the liver lobule in both humans and mice, facilitating a better understanding of the form and regulation of zonation (29).
In the present study, we identified 94 hepatocyte-specified zonation markers (39 central markers and 55 portal markers) by combining four scRNA-seq cohorts for the first time. Based on the zonation markers, we identified three HCC subclasses. Cluster1 was barely involved in zonation-related signature, with a bad prognosis, high recurrence, high AFP level, and advanced TNM stage and histologic grade. Cluster2 was associated with a central signature, with a middle prognosis, AFP level, TNM stage, and histologic grade. Cluster3 exhibited a portal-associated signature, with favorable prognosis, low recurrence, low AFP level, and early TNM stage and histologic grade. In general, this study investigated the zonation characteristics of HCC, and identified three types of HCC with non-zonational-like (Cluster1), central-like (Cluster2), and portal-like (Cluster3) features, respectively.
Moreover, the transcriptome and mutation features, immune infiltration, and immunotherapy response of the subclasses were investigated. Cluster1 (non-zonational-like type) was mainly enriched in differentiation and development-relevant processes, with a high rate of TP53 mutation (42%), a high level of immune cell infiltration, a high expression level of immune checkpoint genes, and may be more sensitive to the CTLA4 inhibitor. Moreover, Cluster2 (central-like type) was mainly enriched in the processes associated with molecular transport and WNT signaling pathway, with a high rate of CTNNB1 mutation (49%), a low level of immune cell infiltration, and a low expression level of immune checkpoint genes. Cluster3 (portal-like type) was enriched in numerous metabolism-associated biological processes and pathways, with low rates of both TP53 (26%) and CTNNB1 (10%), a high level of immune cell infiltration, a high expression level of immune checkpoint genes, and may be more sensitive to the PD1 inhibitor.
Furthermore, we compared the correlation between our clusters and previously published HCC subclasses, which provided us with abundant characteristic information related to our novel zonationassociated clusters, leading to a deeper understanding of our clusters. The results from both RNA-seq and microarray cohorts demonstrated that Cluster1 was associated with the HCC type of Chiang's Interferon class and Deśert's ECM/STEM type. Chiang's Interferon class is featured by the overexpression of several interferon-stimulated genes, with a low mutation rate of CTNNB1 and low expression of IGF2 and CTNNB1 target genes (9). Deśert's ECM/STEM type is the collective name for the ECM type and STEM type. ECM type is characterized by ECM modeling, integrin signaling, and epithelial-mesenchymal transition. STEM type is typified by up-regulation of cell cycle progression and p53 mutation. ECM/STEM-type HCCs share high tumor cell proliferation and a bad prognosis (11). Cluster2 was linked to Deśert's Perivenous type. Deśert's Perivenous type is defined by a high level of perivenous hepatocyte signatures, such as fatty acid and bile salt metabolism, with a high rate of CTNNB1 mutations (11). Cluster3 was associated with Chiang's Proliferation class, Deśert's Periportal type, and Hoshida's S3. Chiang's Proliferation class is characterized by high proliferation, chromosomal instability, activation of IGF and Akt/mTOR signaling, and reduced frequencies of CTNNB1 mutation (9). Deśert's Periportal type is featured by the enrichment of differentiated periportal hepatocyte signatures (gluconeogenesis, amino acid catabolism, hepatocyte nuclear factor 4A (HNF4A) induced genes) and TP53 mutation rates, with good prognosis, low recurrence, early-stage tumors by BCLC, TNM staging systems, and no macrovascular invasion (11). Hoshida's S3 is notable for a well-differentiated hepatocyte signature with a good prognosis, and preserved TP53 function (10). In general, the clusters we constructed validated the findings of previously published HCC subclasses but also preserved their own characteristics. HCC tumors are highly heterogeneous among individuals, and finding targeted therapeutic strategies for specific groups is of vital importance to maximize the therapeutic effect. In the present study, a total of 34 potential drug targets and 16 agents for Cluster1 HCC patients were identified. Among the potential drug targets, AURKB,  (38) in various cancers, including HCC. Therefore, multiple compounds have been designed for these therapeutic targets (Table S11). Among these compounds, Teniposide, K-858, TAK-901, Dexrazoxane, Mitoxantrone, AT-9283, Axitinib, Podophyllotoxin, Amonafide, Amsacrine, Idarubicin, Doxorubicin, Etoposide, GSK1070916, YM−155, and Alisertib were found to be more sensitive for Cluster1 HCC patients.
Teniposide, Dexrazoxane, Mitoxantrone, Podophyllotoxin, Amsacrine, Doxorubicin, and Etoposide are topoisomerase inhibitors that have been widely employed as chemotherapy in multiple cancer therapies. Amonafide is also a topoisomerase II inhibitor and DNA intercalator that has been found to have marked antineoplastic efficacy in preclinical models of cancer. Phase I and II studies revealed that Amonafide might be a promising drug for treating older patients, including those with multidrug-resistant, cytogenetically unfavorable secondary and treatment-associated acute myeloid leukemia (AML) (39).
TAK-901, AT-9283, GSK1070916, and Alisertib are Aurora kinase inhibitors. Human Aurora kinases, including Aurora kinase A (AURKA), B (AURKB), and C (AURKC), play critical roles in monitoring the mitotic checkpoint, bipolar mitotic spindle creation, and centrosome alignment, and also participate in regulating the separation of centrosome, bio-orientation of chromosomes, and cytokinesis. AURKB is a key regulator of mitosis and centrosome through polymerizing microfilaments and regulating chromatid segregation (40). Preclinical studies and Phase I/II/III trials have shown favorable pharmacokinetic properties and markedly antitumor effects of these agents (41)(42)(43)(44).
K-858 is an inhibitor of the mitotic kinesin KIF11 (also known as Eg5). K-858 induces cell mitotic arrest with the formation of monopolar spindles by blocking centrosome separation, and activating the spindle checkpoint (45). K-858 has shown potent anti-proliferative and pro-apoptotic effects in multiple types of cancer cell lines (46,47).
Axitinib is a potent and specific inhibitor of tyrosine kinase, specifically with high affinity for VEGFRs 1, 2, and 3, which was approved by the FDA for use in renal cell carcinoma in 2012 (48).
In addition, Axitinib was also found to inhibit PLK4 with an IC 50 value of 4.2 nM (49). Clinical trials have been performed to treat HCC patients with Axitinib. In a phase II trial, second-line treatment with Axitinib combined with the best supportive care led to a remarkably longer PFS and a higher clinical benefit rate than patients with only best supportive care in advanced HCC (50). Moreover, another multicenter phase II trial in advanced HCC revealed that second-line Axitinib treatment resulted in a 62.2% disease control rate and a 6.7% response rate for advanced HCC after failing the first-line sorafenib therapy (51).
YM-155 is an imidazolium-based compound that has selective activity against BIRC5 (also known as Survivin) (52). Several phase I/II clinical trials have demonstrated that YM-155 has a safe profile and anti-tumor capacity in multiple types of tumors (53,54). Preclinical studies found that YM155 substantially suppressed the proliferation and induced cell cycle arrest and apoptosis of HCC cell lines. Moreover, in a mouse model using patient-derived HCC xenografts with overexpressed BIRC5 and p-BIRC5, YM155 exhibited stronger anti-proliferative efficacy than sorafenib (55).

Conclusion
In conclusion, the current study established a novel HCC classification based on liver zonation signature. By classifying HCC into three clusters with non-zonational-like (Cluster1), central-like (Cluster2), and portal-like (Cluster3) features, this study shed new lights on the heterogeneity of HCC. Cluster3 was intensively correlated with portal-like signature with a good prognosis. Cluster2 was intensively associated with a central-like signature and exhibited low immune infiltration and sensitivity towards immune blockade therapy. Cluster1 exhibited a nonzonational-like signature with the worst prognosis. Moreover, we identified potential drug targets and agents for Cluster1 HCC, which provided new insights into delivering precision medicine for HCC patients and shed new light on improving the therapeutic effect of anti-tumor drugs by selecting potentially sensitive patients.

Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found within the article/Supplementary materials.

Author contributions
TZ, and YN designed this work. JG, XW, and YL integrated and analyzed the data. HL and KC wrote this manuscript. TZ, YN, XC, and JW critical revision of the manuscript for important intellectual content, obtaining funding, and supervision. All authors contributed to the article and approved the submitted version.

Funding
This work was supported by the National Natural Science Foundation of China (NO. 82200972 to TZ, NO. 82070647 to JW, and 82071251 to XC).